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Abstract. We implement a general numerical calculation that allows for a direct comparison between 
nonlinear Hamiltonian dynamics and the Boltzmann-Gibbs canonical distribution in Gibbs -T-space. Us- 
ing paradigmatic first-neighbor models, namely, the inertial XY ferromagnet and the Fermi-Pasta-Ulam 
/9-model, we show that at intermediate energies the Boltzmann-Gibbs equilibrium distribution is a conse- 
quence of Newton second law (F — ma). At higher energies we discuss partial agreement between time 
and ensemble averages. 
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The problem of the dynamical foundation of 
Boltzmann-Gibbs (BG) statistical mechanics dates back 
to the original proposal of this powerful formalism (see, 
PO) despite many important results this fun- 
damental question |5] still presents open basic aspects 
(see, e.g., |HE]E1E1IZ| and references therein). Thanks to 
the current computational capability we can numerically 
integrate the Hamilton equations of large enough sys- 
tems and compare the results with the predictions of the 
BG formalism. Although this technique has been largely 
and successfully implemented in a microcanonical perspec- 
tive (fixed-energy molecular dynamics), the methods used 
when addressing systems in contact with a thermostat 
(such as Monte Carlo and Nose-Hoover [H]) usually impose 
an ad hoc dynamics. In this paper we introduce a scheme 
which enables the discussion of the canonical distribution 
in Gibbs F-space on the basis of the equations of motions. 
Within the present approach both time and ensemble av- 
erages are performed dynamically, so that we are able to 
discuss ergodicity. Using two paradigmatic first-neighbor 
nonlinear Hamiltonian systems — the one-dimensional in- 
ertial XY ferromagnet and the Fermi-Pasta-Ulam (FPU) 
/?-model — we find a remarkable agreement between BG 
equilibrium calculations and dynamical ensemble aver- 
ages. We also compare partial ergodicity failure with the 
maximum Lyapunov coefficient. Our numerical calcula- 
tion can be implemented in systems that allow for a text- 



book definition of canonical ensemble (i.e., part of a large 
isolated system). It would also be interesting to check the 
same procedure in situations where, due for example to 
the presence of long-range terms, important deviations 
from the BG predictions have been found pillO| . We are 
presently making progress on this task. 



Given some macroscopic conditions in the phase space 
of the system under consideration (F-space), the average 
value of a dynamical function can be defined using time or 
ensemble averages; ergodicity means that these two meth- 
ods are equivalent. We remark that both approaches are 
dynamically realizable. In the first case one focuses on a 
single dynamical realization. The probability p/j of find- 
ing the system inside a coarse-grained region R of -T-space 
is defined by the fraction of time tfj spent by the system 
inside that region during the (eventually infinite) total 
amount of time t of its phase space trajectory: = tn/r, 
where the superscript t stands for time definition. The sec- 
ond is achieved for instance by fixing a certain instant of 
time t* and repeating the dynamical evolution up to t*, 
under the same macroscopic (but different microscopic) 
initial conditions. Counting the number of times ur the 
system is found in region R at time t*, with respect to 
the (eventually infinite) total number of times n the cal- 
culation is performed, one defines = nn/n, where the 
superscript e indicates ensemble. 
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For a typical A^-body conservative Hamiltonian sys- 
tem (typical in the sense that it complies with the BG 
prescriptions) at fixed energy (microcanonical setup), 
a standard introduction of the canonical ensemble is ob- 
tained defining the canonical system as composed by a 
subset of M interacting elements, with 1 <C M ^ A^. The 
energy of the M elements satisfies Em ^ E^, and the 
interaction energy between the canonical system and the 
rest of the isolated system (thermal bath) is assumed to 
be much smaller than Em - Under these circumstances, the 
probability pj of finding the system in a M-microstate j 
is given by the BG equilibrium calculation pj cx e~f^^\ 
where (3 = 1/T is the inverse temperature (without loss of 
generality, we set the Boltzmann constant ks = 1), and 
Ej is the energy of the microstate. A dynamical approach 
for the confirmation of this result must face the follow- 
ing numerical difficulty. The /^-space is Md dimensional, 
d being the dimension of the single-particle phase space. 
If we implement a coarse-graining for example by mak- 
ing a partition of k intervals in each coordinate, the total 
number of (hyper)cells is of order k^^'^. Just to put 
some indicative numbers, with fc = 4, A/ = 100 and d = 2 
we get f^M ~ 4^°" ~ 10^^". We should hence implement a 
numerical integration of 2N{':^ 200) Hamilton equations 
with a total amount of time r (or a total number of real- 
izations n) much larger than 10^^*^, which is beyond what 
we can presently do numerically. 

Nevertheless, we can proceed through an alternative 
path and, instead of focusing on the probability associ- 
ated to a microstate, we could consider the probability of 
finding the canonical system with a given energy Em- In 
this case the BG answer is 



io{EM)e 



-13 Em 



p{Em) 

where Z is the partition function and 



(1) 



„ M 

^{Em) = / \{{dpdq,)5[EM - HM{p^,q^)] (2) 
i=l 

is the phase-space density of states at energy Em- As well 
known for a classical system, uj{Em) does not depend on 
any particular thermal statistics, it only depends on the 
Hamiltonian of the system. In other words, we can calcu- 
late u){E]\i) by using any statistics, for example the BG 
one The density of states lu{Em) can be analytically 
estimated through the thermodynamic relation linking the 
statistical entropy to the temperature: d\nuj{E)/dE = p. 
Integrating this relation we have that uj{Em) is given 
through the caloric curve T{E): 



i^{Em) 
uj{Eo) 



exp 



Eo 



dE' (3{E') 



(3) 



where Eq is the energy of the fundamental state. In brief, 
the Hamiltonian structure of the system defines the den- 
sity of states as a function of the energy; once this relation 
is known it is sufficient to multiply uj{Em) by the Boltz- 
mann factor e^P^f-' and to normalize, in order to obtain 



p{Em) for the whole spectrum of temperatures. Now, the 
dynamical computation of p(Em) is much easier than the 
one for pj. All we have to do is to numerically integrate 
Hamilton equations and to calculate the value of the en- 
ergy Em for the canonical subset at each integration step. 
We can then coarse-grain the energy spectrum into bins 
of width AEm and build up a normalized histogram of 
the occurrence of each of these bins. In analogy with the 
previous discussion. 
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represent then the probability distribution of finding the 
canonical system with energy Em , respectively using time 
and ensemble averages. 

To illustrate this calculation, we consider next a spe- 
cific class of analytically solvable nonlinear first-neighbor 
Hamiltonians, 
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(5) 



with periodic boundary conditions (qat+i = qi)- As a first 
case we analyze a one-dimensional chain of rotors with 

y{qi+i - gi) = 1 - cos(g; 



i+l 



), so that the canonical 
coordinates qi G [0, 2-k) and € IR are respectively the 
angular coordinates and the angular momenta of the (unit 
inertial momenta) rotors. This Hamiltonian is an inertial 
version of the classical XY ferromagnetic spin model and 
constitutes a dynamical prototype for spin systems in sta- 
tistical mechanics jHE]- The model is nearly integrable for 
both low and high energies. The former regime is defined 
for T < 0.05 (specific energy e < 0.05) [1] and it is called 
strong coupling regime, for which the rotors constitute a 
set of oscillators almost linearly coupled. The latter is ob- 
tained say for T > 10 (e > 6) 5 , where the rotors are 
almost free (weak coupling regime). For this model, dy- 
namical deviation from BG statistics where detected both 
in the strong and in the weak coupling regimes. Since our 
main scope is to check our calculation scheme in standard 
situations, we will mainly concentrate in the intermediate 
energy range, and discuss partial disagreement that occurs 
at higher energies. The canonical partition function 



» M 

Zm = / Y\_idpidqi)exp[-f3HM{pi,qi)] 



(6) 



gives, for this model, the specific free energy 
/ = -limM^oo[lnZM/(M/3)] (see. e.g., 0): 



f = -T 



ilnr + ln/o(^)+ln27ri 



(7) 



where Io{x) is the modified Bessel function of the 
first kind of order zero. Inversion of the relation 
E{T) = F- TdF/dT furnishes the BG caloric curve r(e), 
where e = limM^oo Em/M. We then rescale the e-axis by 
a factor M and use the fact that the temperature is an 
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Fig. 1. BG analytical canonical prediction for the iner- 
tial XF-ferromagnetic-rotors with M — 100. (a) Loga- 
rithm of the density of states lum{Em)- (b) p{Em) = 
u!m{Em) &qi{—EM /T) / Z , for different temperatures. 



intensive parameter to get T{Em)- The integration in Eq. 

finally gives lo{Em) for any large-but-finite value of M . 
In Fig. ^a) we plot the logarithm of uj{Em) for the first- 
neighbor rotors with M ~ 100 and Fig. ^h) displays BG 
p{Em) for different values of the temperature T (or of the 
specific energy e). We remark that, thanks to the elemen- 
tary properties of the logarithmic function, it is possible 
to implement this calculation for quite large values of Af , 
since one essentially deals with the exponents. 

Because we are interested in very large values of r 
and n, the dynamical integration of Hamilton equations 
has been performed using the 4th order symplectic Neri- 
Yoshida integrator |12| with an iteration parameter that 
assures an energy conservation AEm/Em ~ 10^'^ (a few 
runs with 10~^ showed that 10"'^ is enough for our scopes). 
In particular, we checked that the energy fluctuations of 
the total system introduced by the finite precision of the 
integrator algorithm is order of magnitudes smaller than 
those that one would have in presence of a thermal cou- 
pling. An important point to perform an efficient calcu- 
lation concerns the initial conditions, that must be close 
enough to equilibrium to avoid long transients. In this way 
we focus only on the equilibrium properties of the model, 
ruling out the possible presence of metastable or quasi- 
stationary states appearing with far-from-equilibrium ini- 
tial conditions. Since the system does not display any 
phase transition for T > but presents a tendency to clus- 
terization at low temperatures, we have used a Maxwellian 
distribution for the angular momenta (with the appro- 
priate temperature) and a set of I equidistant Gaussian 
distributions for the angles, each with the same variance 
appropriately calculated in order to yield the desired to- 
tal energy Eiq. In our calculation it was sufficient to use 
I — 6 for a fast enough relaxation to equilibrium in all 
our (microcanonical) setups. We have checked that this 
particular choice has no influence on the functional form 
of the equilibrium probability density functions: it is done 
to save computational time. Different close-to-equilibrium 
initial conditions eventually yield the same results. For all 
our results we have waited for 10'^ iteration steps before 
starting the measurements for a canonical system which 
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Fig. 2. (a-c) Comparison between the BG prediction p{Em) 
(full line), the ensemble dynamical average p''{Em) (crosses), 
and the time dynamical average p ^(Em) (circles), k = Km /M 
is the value of the specific kinetic energy, (d) Analysis of dis- 
crepancy between p "^(Sm) and p(iJA/) (crosses), andp*(_EM) 
and p{Em) (empty circles). We also plot the largest Lyapunov 
coefficient Xmax (squares) and the inverse of the time-scale for 
a normal fluctuation I/Iaem (full circles). The lines are guides 
to the eye. See text for details. 



is composed by a randomly chosen subset of M adjacent 
rotors. 

In Fig. [21 a-c) we present a striking agreement between 
the BG analytical prediction for p{Em) (full line) and the 
dynamical estimation oi p''{Em) (crosses) for various or- 
der of magnitudes of the specific energy e with a setup 
{M,N) = (10^,10^) and a total number of realizations 
n = 5 X 10^. On the other hand, p*{Em) (circles), calcu- 
lated with a total number of iteration steps r = 5 x 10^, 
displays a good agreement with respect to the BG an- 
alytical distribution for intermediate energies, but starts 
to show large discrepancies when entering in the weak- 
coupling regime. In order to quantify this difference, we 
have defined the discrepancy < e < 2 between two distri- 
butions as the integral of the absolute value of the differ- 
ence of the distributions. To allow for a comparison with 
the largest Lyapunov coefficient, in Fig. |3d) we plot the 
quantity < 1/e — 1/2 < 00 which is zero for maximum 
discrepancy and infinite for perfect overlap of the distri- 
butions. While for ensemble averages 1/e — 1/2 is large 
and almost constant with the energy, in the case of time 
averages such a quantity presents a dramatic decrease for 
large energies. In fact, we verified that the time necessary 
to have a typical energy fluctuation of the canonical subset 
(AEm ~ Em/VM) grows with the energy (see full circles 
in Fig. Efd) , where we plot the inverse of this time) , as 
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Fig. 3. Dynamical evidence of the Boltzmann factor. We plot 
ln[p '^(i?M)/w(-EA/)] for the ensemble averages of Fig. |5| (cir- 
cles). T is the reciprocal of the slope of linear regressions (full 
lines) on the data. Insets (a) and (b) show a magnification of 
the results for e — 0.05 and e = 0.5 respectively. 



a consequence of the fact that rotors are increasingly free 
(the potential is upper bounded). We point out that the 
largest Lyapunov coefficient (squares in Fig.l^d)) does not 
display a significant correlation with the time character- 
izing relaxation of p^{Em) (circles in Fig.|2Ia-c)) towards 
the BG distribution p{Em) (see also |5| for a discussion 
of this point). This means that for the present system the 
positivity of the largest Lyapunov coefficient is a measure 
of local chaos that does not imply relaxation to global 
chaos. 

An important result is the coincidence between the 
value of the BG temperature T and twice the specific ki- 
netic energy k = Km/M within an error of at most 2%. 
We stress that the probability density functions shown in 
Fig. 121 are obtained by means of first principles only and 
with complete independence from the BG theory, which 
we are checking. The concurrence between dynamics and 
the Boltzmann factor appears satisfactorily in the linear 
regressions of Fig. 13 where we plot ln[p '^(£'M)/t^(£'A/)] 
for the ensemble averages of Fig. Ela-c). With other val- 
ues of (M, iV), namely (50, 500) and (10^, lO^), the results 
were qualitatively the same. 

We also obtained a confirmation of our results by im- 
plementing the same calculation scheme for the FPU (3- 
model, defined by the potential V{qi+x — qt) = {qi+i — 
qif' jl -t- 0.1(gi+i — (?i)*/4 with qi G R, again considering 
close-to-equilibrium initial conditions (see, e.g., for the 
analytical canonical solution and for a discussion of initial 
conditions). Although it is known that the FPU model 
presents, in common with the rotors model, a very rich 
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Fig. 4. Same as Fig. |2a-c) for the FPU /3- model. See text for 
details. 



anomalous behavior at low energies for our initial 

conditions and for the energy-range we tested we found 
that p*{Em) is in good agreement with the BG predic- 
tion (Fig. Ell- 
in summary, we recall that using the standard BG for- 
malism and common numerical techniques, we have in- 
troduced a new calculation that allows for a comparison 
between nonlinear Newtonian dynamics and canonical sta- 
tistical mechanics. Implementing a standard setup we have 
in fact shown that the BG energy distribution in .T-space 
coincides with the one that is obtained dynamically (in- 
tegrating Hamilton equations for close-to-equilibrium ini- 
tial conditions) when an ensemble average is executed. We 
have checked this conclusion for two paradigmatic first- 
neighbor nonlinear Hamiltonians. As a side result, this 
calculation provides a dynamical confirmation of the very 
well known relation between temperature and specific ki- 
netic energy k — T/2 (for one-dimensional systems). With 
respect to finite-time dynamical averages, at moderate low 
energies we have found a confirmation of the BG predic- 
tions. For the XY-mode\ at high energies, if the time- 
scale is not very large, finite-time averages disagree with 
ensemble averages as a consequence of an increase of the 
time-scale of typical energy fluctuations. The energy de- 
pendence of this discrepancy does not display correlation 
with that of the largest Lyapunov coefficient (see also [5]). 

Finally, let us emphasize that what we have shown 
here is that equilibrium thermal statistics descends from 
(finite-precision) mechanics, even for a system in contact 
with a thermostat (usually discussed through Monte Carlo 
or Nosc-Hoover techniques, which do not deduce the equi- 
librium distribution but impose it |H]). Indeed, this is the 
significance of Figs.[2Ia-c) andQ] where circles and crosses 
have been obtained from Newton law, whereas full lines 
come from the BG theory. Equivalently, if we recall that 
the density of states is a purely mechanical concept, the 
same conclusion is shown in Fig. 13 The present calcula- 
tion scheme provides an insight onto the basic question of 
the dynamical foundation of statistical mechanics ^|21|51 

J and may serve as a useful tool in the discussion of 
complex situations (see e.g., (Hj) where dynamical discrep- 
ancies with the BG theory have been found. 
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